In [131]:
import sklearn
import sklearn.datasets
import sklearn.ensemble
import numpy as np
import pandas as pd
import lime
import lime.lime_tabular
from __future__ import print_function
np.random.seed(1)

Numerical and Categorical features in the same dataset¶

We will analyse a dataset that has both numerical and categorical features. Here, the task is to predict whether a person makes over 50K dollars per year.

In [132]:
feature_names = ["Age", "Workclass", "fnlwgt", "Education", "Education-Num", "Marital Status", "Occupation", "Relationship", "Race", "Sex", 
                 "Capital Gain", "Capital Loss","Hours per week", "Country"]
In [133]:
data = np.genfromtxt('https://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data', delimiter=', ', dtype=str)

Take a look at the data. Let's analyse the labels:

In [134]:
labels = data[:, 14]
np.unique(labels)
Out[134]:
array(['<=50K', '>50K'], dtype='<U26')

In order to use it, we need to preprocess the labels to have discrete values:

In [135]:
le=sklearn.preprocessing.LabelEncoder()
le.fit(labels)
labels = le.transform(labels)
class_names = le.classes_
In [136]:
np.unique(labels)
Out[136]:
array([0, 1])

Let's remove labels from data and take a look at the data:

In [137]:
data = data[:, :-1]
In [138]:
pd.DataFrame(data, columns=feature_names)
Out[138]:
Age Workclass fnlwgt Education Education-Num Marital Status Occupation Relationship Race Sex Capital Gain Capital Loss Hours per week Country
0 39 State-gov 77516 Bachelors 13 Never-married Adm-clerical Not-in-family White Male 2174 0 40 United-States
1 50 Self-emp-not-inc 83311 Bachelors 13 Married-civ-spouse Exec-managerial Husband White Male 0 0 13 United-States
2 38 Private 215646 HS-grad 9 Divorced Handlers-cleaners Not-in-family White Male 0 0 40 United-States
3 53 Private 234721 11th 7 Married-civ-spouse Handlers-cleaners Husband Black Male 0 0 40 United-States
4 28 Private 338409 Bachelors 13 Married-civ-spouse Prof-specialty Wife Black Female 0 0 40 Cuba
... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
32556 27 Private 257302 Assoc-acdm 12 Married-civ-spouse Tech-support Wife White Female 0 0 38 United-States
32557 40 Private 154374 HS-grad 9 Married-civ-spouse Machine-op-inspct Husband White Male 0 0 40 United-States
32558 58 Private 151910 HS-grad 9 Widowed Adm-clerical Unmarried White Female 0 0 40 United-States
32559 22 Private 201490 HS-grad 9 Never-married Adm-clerical Own-child White Male 0 0 20 United-States
32560 52 Self-emp-inc 287927 HS-grad 9 Married-civ-spouse Exec-managerial Wife White Female 15024 0 40 United-States

32561 rows × 14 columns

The dataset has many categorical features, which we need to preprocess like we did with the labels before - our explainer (and most classifiers) takes in numerical data, even if the features are categorical. We thus transform all of the string attributes into integers, using sklearn's LabelEncoder. We use a dict to save the correspondence between the integer values and the original strings, so that we can present this later in the explanations.

In [139]:
categorical_features = [1, 3, 5, 6, 7, 8, 9, 13]
In [140]:
categorical_names = {}
for feature in categorical_features:
    le = sklearn.preprocessing.LabelEncoder()
    le.fit(data[:, feature])
    data[:, feature] = le.transform(data[:, feature])
    categorical_names[feature] = le.classes_
In [141]:
data = data.astype(float)

Final look at the preprocessed data:

In [142]:
pd.DataFrame(data, columns=feature_names)
Out[142]:
Age Workclass fnlwgt Education Education-Num Marital Status Occupation Relationship Race Sex Capital Gain Capital Loss Hours per week Country
0 39.0 7.0 77516.0 9.0 13.0 4.0 1.0 1.0 4.0 1.0 2174.0 0.0 40.0 39.0
1 50.0 6.0 83311.0 9.0 13.0 2.0 4.0 0.0 4.0 1.0 0.0 0.0 13.0 39.0
2 38.0 4.0 215646.0 11.0 9.0 0.0 6.0 1.0 4.0 1.0 0.0 0.0 40.0 39.0
3 53.0 4.0 234721.0 1.0 7.0 2.0 6.0 0.0 2.0 1.0 0.0 0.0 40.0 39.0
4 28.0 4.0 338409.0 9.0 13.0 2.0 10.0 5.0 2.0 0.0 0.0 0.0 40.0 5.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
32556 27.0 4.0 257302.0 7.0 12.0 2.0 13.0 5.0 4.0 0.0 0.0 0.0 38.0 39.0
32557 40.0 4.0 154374.0 11.0 9.0 2.0 7.0 0.0 4.0 1.0 0.0 0.0 40.0 39.0
32558 58.0 4.0 151910.0 11.0 9.0 6.0 1.0 4.0 4.0 0.0 0.0 0.0 40.0 39.0
32559 22.0 4.0 201490.0 11.0 9.0 4.0 1.0 3.0 4.0 1.0 0.0 0.0 20.0 39.0
32560 52.0 5.0 287927.0 11.0 9.0 2.0 4.0 5.0 4.0 0.0 15024.0 0.0 40.0 39.0

32561 rows × 14 columns

In [143]:
pd.DataFrame(data, columns=feature_names).describe()
Out[143]:
Age Workclass fnlwgt Education Education-Num Marital Status Occupation Relationship Race Sex Capital Gain Capital Loss Hours per week Country
count 32561.000000 32561.000000 3.256100e+04 32561.000000 32561.000000 32561.000000 32561.000000 32561.000000 32561.000000 32561.000000 32561.000000 32561.000000 32561.000000 32561.000000
mean 38.581647 3.868892 1.897784e+05 10.298210 10.080679 2.611836 6.572740 1.446362 3.665858 0.669205 1077.648844 87.303830 40.437456 36.718866
std 13.640433 1.455960 1.055500e+05 3.870264 2.572720 1.506222 4.228857 1.606771 0.848806 0.470506 7385.292085 402.960219 12.347429 7.823782
min 17.000000 0.000000 1.228500e+04 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000
25% 28.000000 4.000000 1.178270e+05 9.000000 9.000000 2.000000 3.000000 0.000000 4.000000 0.000000 0.000000 0.000000 40.000000 39.000000
50% 37.000000 4.000000 1.783560e+05 11.000000 10.000000 2.000000 7.000000 1.000000 4.000000 1.000000 0.000000 0.000000 40.000000 39.000000
75% 48.000000 4.000000 2.370510e+05 12.000000 12.000000 4.000000 10.000000 3.000000 4.000000 1.000000 0.000000 0.000000 45.000000 39.000000
max 90.000000 8.000000 1.484705e+06 15.000000 16.000000 6.000000 14.000000 5.000000 4.000000 1.000000 99999.000000 4356.000000 99.000000 41.000000

As we see, the categorical data has numerical values indicating the categories now.
We now split the data into training and testing:

In [144]:
train, test, labels_train, labels_test = sklearn.model_selection.train_test_split(data, labels, train_size=0.80)

Finally, we use a One-hot encoder, so that the classifier does not take our categorical features as continuous features. We will use this encoder only for the classifier, not for the explainer - and the reason is that the explainer must make sure that a categorical feature only has one value set to True.

In [145]:
from sklearn.compose import ColumnTransformer 
from sklearn.preprocessing import OneHotEncoder
In [146]:
encoder = ColumnTransformer([("OneHot", OneHotEncoder(), categorical_features)], remainder = 'passthrough')
In [147]:
encoder.fit(data)
encoded_train = encoder.transform(train)

We will use gradient boosted trees as the model, using the xgboost package.

In [148]:
import xgboost
gbtree = xgboost.XGBClassifier(n_estimators=300, max_depth=5)
gbtree.fit(encoded_train, labels_train)
Out[148]:
XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, device=None, early_stopping_rounds=None,
              enable_categorical=False, eval_metric=None, feature_types=None,
              gamma=None, grow_policy=None, importance_type=None,
              interaction_constraints=None, learning_rate=None, max_bin=None,
              max_cat_threshold=None, max_cat_to_onehot=None,
              max_delta_step=None, max_depth=5, max_leaves=None,
              min_child_weight=None, missing=nan, monotone_constraints=None,
              multi_strategy=None, n_estimators=300, n_jobs=None,
              num_parallel_tree=None, random_state=None, ...)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, device=None, early_stopping_rounds=None,
              enable_categorical=False, eval_metric=None, feature_types=None,
              gamma=None, grow_policy=None, importance_type=None,
              interaction_constraints=None, learning_rate=None, max_bin=None,
              max_cat_threshold=None, max_cat_to_onehot=None,
              max_delta_step=None, max_depth=5, max_leaves=None,
              min_child_weight=None, missing=nan, monotone_constraints=None,
              multi_strategy=None, n_estimators=300, n_jobs=None,
              num_parallel_tree=None, random_state=None, ...)
In [149]:
sklearn.metrics.accuracy_score(labels_test, gbtree.predict(encoder.transform(test)))
Out[149]:
0.868417012129587

Our predict function, which first transforms the data into the one-hot representation:

In [150]:
predict_fn = lambda x: gbtree.predict_proba(encoder.transform(x)).astype(float)

Explaining predictions¶

Tabular explainers need a training set. The reason for this is because we compute statistics on each feature (column). If the feature is numerical, we compute the mean and std, and discretize it into quartiles. If the feature is categorical, we compute the frequency of each value. For this tutorial, we'll only look at numerical features.

We use these computed statistics for two things:

  • To scale the data, so that we can meaningfully compute distances when the attributes are not on the same scale
  • To sample perturbed instances - which we do by sampling from a Normal(0,1), multiplying by the std and adding back the mean.

We now create our explainer. The categorical_features parameter lets it know which features are categorical. The categorical names parameter gives a string representation of each categorical feature's numerical value.

In [151]:
explainer = lime.lime_tabular.LimeTabularExplainer(train, feature_names=feature_names, class_names=class_names,
                                                   categorical_features=categorical_features, 
                                                   categorical_names=categorical_names, kernel_width=3, verbose=True)

We now show a few explanations with a verbose set to True.

In [152]:
np.random.seed(1)
i = 1653
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('xgboost_output.html')
Intercept -0.004040331593255592
Prediction_local [1.03542247]
Right: 0.9999810457229614
In [153]:
exp.as_list()
Out[153]:
[('Capital Gain > 0.00', 0.7137120588369433),
 ('Marital Status=Married-civ-spouse', 0.10496853816074189),
 ('Education-Num > 12.00', 0.08327993642177214),
 ('Hours per week > 45.00', 0.07942690425041142),
 ('Age > 48.00', 0.05807536707034083)]

First, note that the row we explained is displayed on the right side, in table format. Since we had the show_all parameter set to false, only the features used in the explanation are displayed.

The "value" column displays the original value for each feature.

The explanations for categorical features are based not only on features, but on feature-value pairs.

LIME has discretized the features in the explanation. This is because we let discretize_continuous=True in the constructor (this is the default). Discretized features make for more intuitive explanations.

As for the values displayed after setting the "verbose" parameter to True: Intercept is the intercept of the linear model used inside the LIME algorithm. Prediction_local is the prediction of this model for the instance of interest, and Right is the xgboost model's prediction for the same instance. We analyse the weights with respect to the intercept.

Note that capital gain has very high weight. This makes sense. Now let's see an example where the prediction is different:

In [154]:
i = 92
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('xgboost_output2.html')
Intercept 0.7999674614920426
Prediction_local [0.10059035]
Right: 0.07562913000583649

Let's also analyse one "weird" example. Take a look at the explanations:

In [155]:
i = 18
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('xgboost_output3.html')
Intercept 0.12810762181146226
Prediction_local [0.83525275]
Right: 0.005185970105230808

We see that the model predicted the output "<=50K" even though the explanation suggests a different result. Let's do the analysis:
intercept + weights = prediction_local, which is our linear model's prediction, corresponding to label ">50K". The xgboost's prediction is, however, 0.005, so in this case the LIME explainer is not useful.

Homework¶

  • Choose a different model for preparing predictions (you may want to use sklearn models).
  • Prepare an explainer.
  • Use it on the three instances explained in the tutorial. Did the explanations change? Analyse if the explanations make sense.
  • Explain an example where the model prediction was incorrect. Comment on the results.
  • Analyse one more example with an interesting explanation.
In [156]:
from sklearn.ensemble import AdaBoostClassifier
from sklearn.tree import DecisionTreeClassifier

base_estimator = DecisionTreeClassifier(max_depth=1)

ada = AdaBoostClassifier(estimator=base_estimator, n_estimators=50, learning_rate=1.0, random_state=42, algorithm='SAMME')

ada.fit(encoded_train, labels_train)
Out[156]:
AdaBoostClassifier(algorithm='SAMME',
                   estimator=DecisionTreeClassifier(max_depth=1),
                   random_state=42)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
AdaBoostClassifier(algorithm='SAMME',
                   estimator=DecisionTreeClassifier(max_depth=1),
                   random_state=42)
DecisionTreeClassifier(max_depth=1)
DecisionTreeClassifier(max_depth=1)
In [157]:
print(f'accuracy {sklearn.metrics.accuracy_score(labels_test, ada.predict(encoder.transform(test)))}')
accuracy 0.857976354982343
In [158]:
predict_fn = lambda x: ada.predict_proba(encoder.transform(x)).astype(float)
In [159]:
explainer = lime.lime_tabular.LimeTabularExplainer(train, feature_names=feature_names, class_names=class_names,
                                                   categorical_features=categorical_features, 
                                                   categorical_names=categorical_names, kernel_width=3, verbose=True)

First example >50¶

In [178]:
np.random.seed(3)
i = 1653
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('ada_output1.html')
Intercept 0.2878457554605127
Prediction_local [0.72477611]
Right: 0.7573225370099033

Second example <=50¶

In [163]:
i = 92
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('ada_output2.html')
Intercept 0.5935723418545474
Prediction_local [0.37835424]
Right: 0.43384415942339394

Third example - weird¶

In [190]:
i = 18
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('ada_output3.html')
Intercept 0.38669199559782574
Prediction_local [0.56816298]
Right: 0.43081361135082946

Explanation for AdaBoost are another than for XGBoost. Accuracy is insignificantly lower, but prediction probability for given examples are significantly lower. Generally the same features (in xgboost and ada) were chosen as important to predict these examples.

Fourth example - incorrect classification, predicted <=50k¶

In [206]:
i = 5
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('ada_output4.html')
Intercept 0.5084189740754964
Prediction_local [0.4741731]
Right: 0.48551043153409934
In [205]:
i = 5
print('true: ',  labels_test[i])
print('predicted: ', [0 if predict_fn(test[i].reshape(1, -1))[0][0]>0.5 else 1][0])
true:  1
predicted:  0

In this example model has classified incorrectly, but had a small assurance in it's decision. It classified class <=50k with 51% probability. Small capital gain was argued to class <=50k even more than fact he is married and working over 45h/week was argued to class >50k.

Fifth example - incorrect classification, predicted >50k¶

In [219]:
for i in range (1, len(labels_test)):
    if labels_test[i] == 0 and [0 if predict_fn(test[i].reshape(1, -1))[0][0]>0.5 else 1][0] == 1:
        print(i)
        break
43
In [220]:
i = 43
print('true: ',  labels_test[i])
print('predicted: ', [0 if predict_fn(test[i].reshape(1, -1))[0][0]>0.5 else 1][0])
true:  0
predicted:  1
In [221]:
i = 43
exp = explainer.explain_instance(test[i], predict_fn, num_features=5)
exp.show_in_notebook(show_all=False)
exp.save_to_file('ada_output5.html')
Intercept 0.5233645381369108
Prediction_local [0.49205317]
Right: 0.5102165703703702

Again model wasn't sure in decision. The correct class was <=50k and model has predicted >50k with 51% probability. Small capital gain and capital loss couldn't argued to <=50k more than fact about marriage and working less than 45h/week.

Comparing 4 and 5 examples it can be observed that tiny changes in data can change predicted class. The top 5 features didn't decide about predicted class.

In [ ]: